beagle<-read.table("beagle_chr20_util_gen_r2",sep="\t")
impute<-read.table("impute_chr20_util_gen_r2",sep="\t")
colnames(beagle)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
colnames(impute)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
merged<-merge(beagle,impute,by="marker_identifier")
dim(merged)
beagle1<-read.table("beagle_out_20.snps_chr20.bgl.gz.r2",sep="\t")
dim(beagle1)
colnames(beagle1)<-c("marker_identifier","r2")
merged<-merge(merged,beagle1,by="marker_identifier")
dim(merged)
merged_f<-merged[which(merged$r2 !='NaN'),]
dim(merged_f)
merged_f[1,]
main<-paste("Beagle allelic r2 VS Impute allelic r2,pearson corr=",cor(merged_f$allelic_r_squared.x,merged_f$allelic_r_squared.y,method="pearson"),sep=" ")
main
pdf("beagle_allelic_r2_vs_impute_allelic_r2_chr20.pdf")
plot(merged_f$allelic_r_squared.x,merged_f$allelic_r_squared.y,xlab = c("beagle allelic r2"),ylab=c("impute allelic r2"),main=main,pch=".")
dev.off()
impute<-read.table("chr20.one.phased.impute2_fil_chr20_info",sep=" ",head=T)
dim(impute)
colnames(impute)<-c("marker_identifier","position")
merged_f<-merge(merged_f,impute,by="marker_identifier")
dim(merged_f)
impute[1,]
merged_f<-merged[which(merged$r2 !='NaN'),]
dim(merged_f)
impute<-read.table("chr8.one.phased.impute2_fil_chr8_info",sep=" ",head=T)
impuete[1,]
impute[1,]
impute<-impute[,c(2,3)]
merged_f[1,]
colnames(impute)<-c("marker_identifier","position")
merged_f<-merge(merged_f,impute,by="marker_identifier")
dim(merged_f)
dim(impute)
merged_f[1,]
pdf("beagle_allelic_R2VSimpute_allelic_R2.pdf")
plot(merged_f$position,merged_f$allelic_r_squared.x,ylim=c(0,1),xlab = c("POSITION"),ylab = c("Allelic R2 value"),main=c("beagle allelic r2 (red) vs impute allelic r2 (green)") ,pch=".",col=c("red"))
points(merged_f$position,merged_f$allelic_r_squared.y,pch=".",col=c("green"))
dev.off()
multhist <- function(x,main,x.loc,...) {
  new.x <- lapply(x,hist,plot=F)
  tmp <- lapply(new.x,function(x) x$breaks)
  max.break <- lapply(tmp,function(x) max(x))
  max.break <- which(unlist(max.break)==max(unlist(max.break)))[1]
  breaks <- new.x[[max.break]]$breaks
  new.x <- lapply(x,hist,plot=F,breaks=breaks)
  y.max <- max(unlist(lapply(new.x,function(x) max(x$counts))))
  barplot(do.call(rbind,lapply(new.x,function(x) x$counts)),beside=T,
          space=c(0,0),names.arg=signif(new.x[[1]]$mids,2),
          axes=T,legend.text=T,main=main,args.legend=list(x=x.loc,y=y.max),...)
}
pdf("Multiple_hist_chr8_allelic_r2.pdf")
multhist(x=list(Beagle_allelic_r2= merged_f[,4],Impute_allelic_r2= merged_f[,11]),main="Beagle_allelic_r2 VS Impute_allelic_r2", x.loc=20)
dev.off()
savehistory("plot_chr20_allelic.R")
